function [M2, coefs, varargout] = rescaleMatrix(M, varargin)
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% This function takes a matrix as input where each column is a variable, and each
	% column is an observation.
	% It rescales the variables to avoid very large numbers or very small numbers.
	% To do that, it multiplies or divides the whole column by a power of 10.
	% That power of 10 is obtained by taking the max absolute value of the variable,
	%	then take round(log10(x)).
	%
	% Note: M2 = M  .* repmat(coefs, NumRows, 1);
	%  or   M =  M2 ./ repmat(coefs, NumRows, 1);
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%% Inputs:
	% M:			NumRows x NumCols
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%% Outputs:
	% M2:			NumRows x NumCols
	% coefs:		1 x NumCols
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
	% Read dimensions
	[NumRows, NumCols] = size(M);
	
	% Figure out "best" coefs
	tmp = round(log(max(abs(M)))/log(10));	% 1 x NumCols
	coefs = 10.^(-tmp);						% 1 x NumCols
	
	% Get rescaled matrix
	M2 = M .* repmat(coefs, NumRows, 1);	% NumRows x NumCols
	
	% Get new variable names
	if nargin >= 2
		VarNames2 = varargin{1};
		for cc = 1:NumCols
			if coefs(cc) ~= 1
				VarNames2{cc} = sprintf('(%s x 1e%d)', VarNames2{cc}, -tmp(cc));
			end
		end
		varargout{1} = VarNames2;
	end
end
